Travel Time Matrix to Maps

Notebook Purpose

This notebook serves to summarize the entire visualization process going from the travel time matrix to the visualizations. That involves the following sections:

1. Importing and Linking Geospatial Data
2. Interactive Visualization
3. Map HTML Exports

0) Useful Libraries

# import custom scoring, cleaning, and visualization functions
source('../0 - custom_functions/functions.R')

# wrangling/convenience
library(tidyverse)
library(glue)
library(stringr)
library(sf)
library(data.table)
library(bit64)

# visualization
library(leaflet)
library(mapview); mapviewOptions(platform = 'leafgl')

# For pretty knitting
library(lemon)
knit_print.data.frame <- lemon_print
knit_print.tbl <- lemon_print
knit_print.summary <- lemon_print

shapepath <- "../../data/1_raw/shape_files"
comppath <- "../../data/3_computed/"
mappath <- "../../data/html_maps"    # output directory

1) Importing Data

Import the Accessibility Measures

scores_frame <- read.csv(paste0(comppath, '/accessibility_measures/scores_frame.csv'))
isochr_frame <- read.csv(paste0(comppath, '/accessibility_measures/isochrone_frame.csv'))

# convert factor columns to factor
scores_frame[, c(1,2,4,5)] <- lapply(scores_frame[, c(1,2,4,5)], as.factor)
isochr_frame[, c(1,2,3)] <- lapply(isochr_frame[, c(1,2,3)], as.factor)

head(scores_frame)

head(isochr_frame)

NA

Import the dissemination block shape file

Note: The shape file being imported was preproccessed from a national census shape file and does not need to be filtered/cleaned any further. GitBash scripting (via jq library) was used to extract 36 megabytes of Vancouver polygon data from 1.6 gigabytes of Canada wide geoJson polygon data, which R nor Python could efficiently handle.

vancouver_shape <- st_read(paste0(shapepath, '/DB_Van_CMA.shp'), stringsAsFactors = FALSE)

Reading layer `DB_Van_CMA' from data source `C:\Users\Luka\Documents\GitHub\Capstone-Project\data\1_raw\shape_files\DB_Van_CMA.shp' using driver `ESRI Shapefile'
Simple feature collection with 15197 features and 27 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 4001643 ymin: 1957237 xmax: 4068119 ymax: 2032663
Projected CRS: PCS_Lambert_Conformal_Conic

# id to factor
vancouver_shape$DBUID <- as.factor(vancouver_shape$DBUID)
head(vancouver_shape)

Simple feature collection with 6 features and 27 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 4019407 ymin: 2002663 xmax: 4033248 ymax: 2007975
Projected CRS: PCS_Lambert_Conformal_Conic
        DBUID DBRPLAMX DBRPLAMY PRUID                                  PRNAME CDUID            CDNAME CDTYPE  CCSUID
1 59150244005  4030984  2005717    59 British Columbia / Colombie-Britannique  5915 Greater Vancouver     RD 5915020
2 59150244008  4032066  2007480    59 British Columbia / Colombie-Britannique  5915 Greater Vancouver     RD 5915020
3 59150372006  4019917  2002902    59 British Columbia / Colombie-Britannique  5915 Greater Vancouver     RD 5915022
4 59150372007  4019852  2002834    59 British Columbia / Colombie-Britannique  5915 Greater Vancouver     RD 5915022
5 59150373001  4019547  2002835    59 British Columbia / Colombie-Britannique  5915 Greater Vancouver     RD 5915022
6 59150373002  4019483  2002723    59 British Columbia / Colombie-Britannique  5915 Greater Vancouver     RD 5915022
              CCSNAME  CSDUID         CSDNAME CSDTYPE ERUID                                                ERNAME FEDUID
1 Greater Vancouver A 5915046 North Vancouver      DM  5920 Lower Mainland--Southwest / Lower Mainland--Sud-ouest  59002
2 Greater Vancouver A 5915046 North Vancouver      DM  5920 Lower Mainland--Southwest / Lower Mainland--Sud-ouest  59002
3           Vancouver 5915022       Vancouver      CY  5920 Lower Mainland--Southwest / Lower Mainland--Sud-ouest  59038
4           Vancouver 5915022       Vancouver      CY  5920 Lower Mainland--Southwest / Lower Mainland--Sud-ouest  59038
5           Vancouver 5915022       Vancouver      CY  5920 Lower Mainland--Southwest / Lower Mainland--Sud-ouest  59038
6           Vancouver 5915022       Vancouver      CY  5920 Lower Mainland--Southwest / Lower Mainland--Sud-ouest  59038
                                         FEDNAME SACCODE SACTYPE CMAUID CMAPUID   CMANAME CMATYPE      CTUID  CTNAME
1 Burnaby North--Seymour / Burnaby-Nord--Seymour     933       1    933   59933 Vancouver       B 9330110.03 0110.03
2 Burnaby North--Seymour / Burnaby-Nord--Seymour     933       1    933   59933 Vancouver       B 9330110.03 0110.03
3                             Vancouver Kingsway     933       1    933   59933 Vancouver       B 9330032.00 0032.00
4                             Vancouver Kingsway     933       1    933   59933 Vancouver       B 9330032.00 0032.00
5                             Vancouver Kingsway     933       1    933   59933 Vancouver       B 9330032.00 0032.00
6                             Vancouver Kingsway     933       1    933   59933 Vancouver       B 9330032.00 0032.00
    ADAUID    DAUID                       geometry
1 59150014 59150244 MULTIPOLYGON (((4031013 200...
2 59150014 59150244 MULTIPOLYGON (((4031720 200...
3 59150082 59150372 MULTIPOLYGON (((4019999 200...
4 59150082 59150372 MULTIPOLYGON (((4019845 200...
5 59150082 59150373 MULTIPOLYGON (((4019658 200...
6 59150082 59150373 MULTIPOLYGON (((4019485 200...

Import bus stop network data

# Read Vancouver Station txt file
rawlocs <- read.csv("../../data/1_raw/transit_and_osm_data/stops.txt",
                    head = TRUE, sep=",")

station_locs <- rawlocs[c("stop_id","stop_name","stop_lat","stop_lon")]

colnames(station_locs)[3] <- "latitude"
colnames(station_locs)[4] <- "longitude"

# Convert the columns imported as a factor to characters
station_locs$stop_id <- as.character(station_locs$stop_id)
station_locs$stop_name <- as.character(station_locs$stop_name)

# check for NA rows
station_locs[is.na(as.numeric(station_locs$stop_id)), ]

# check how many stops there are
uniqueN(station_locs$stop_id)

[1] 8841

head(station_locs)

NA

Joining the GeoSpatial Data

# merges shape and visualization data, then transforms them for visualizations
# crs = target coordinate reference system
mapping_data_prepper <- function(shape_data, visualization_data,
                                 by = NULL, crs = 4326) {
  if (is.null(by)) {
    message("Must provide key or columns to join on. For example, on = c('ID' = 'id'))")
    return()
  }
  
  # join both datasets on the shapes so all polygons are included
  shape_viz_frame <- left_join(shape_data, visualization_data, by = by)
  
  # convert to simple feature object (sf) then transform coords
  return(st_transform(st_as_sf(shape_viz_frame), crs = crs))
    
}

scores_viz_frame <- mapping_data_prepper(vancouver_shape, scores_frame,
                                         by = c('DBUID' = 'fromId'))

isochr_viz_frame <- mapping_data_prepper(vancouver_shape, isochr_frame,
                                         by = c('DBUID' = 'fromId'))

efficiency_frame <- read.csv('../../data/3_computed/transit_efficiency/efficiency_frame.csv')
efficiency_frame$fromId <- as.factor(efficiency_frame$fromId)
efficiency_frame_viz <- mapping_data_prepper(vancouver_shape, efficiency_frame,
                                             by = c('DBUID' = 'fromId'))

2) Interactive Visualization

# this cell used to be used for visualization experimentation
# the functions have since been migrated to the visualization section of:
# 0 - custom_functions

3) Map HTML Exports

Maps to export: 1) Score maps (32) 2) Isochrone maps (4) 3) Kepler maps (4) - performed elsewhere 4) Efficiency maps (2)

# 2 efficiency maps
map_maker_efficiency_cont(efficiency_frame_viz, 
                          bus_data = station_locs,
                          view = TRUE,
                          output_dir = paste0(mappath, '/efficiency_maps'))

[1] "Current Map: Continuous Efficiency Map"
[1] "Current Map: Continuous Efficiency Map with bus stops"

map_maker_efficiency_cont(efficiency_frame_viz, output_dir = paste0(mappath, '/efficiency_maps'))